import ekf_testv6 as ekf
import numpy as np
# import pykalmani
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import scipy.stats as ss
import matplotlib.transforms as mtransforms
plt.rcParams["font.family"] = "Arial"

def cdf_subplot(ax,loc,lims,xs,ps,pcolor='purple'):
    #loc: threshold, thickness, start (years), width (in years)
    delta=20
    ax.plot([lims[2],loc[2]-delta], [loc[0], loc[0]], '-',color='sienna')
    ax.plot([loc[2]-delta,loc[2]], [loc[0], loc[0]+loc[1]/2.0], '-', color='lightcoral')
    ax.plot([loc[2]-delta,loc[2]], [loc[0], loc[0]-loc[1]/2.0], '-', color='olive')

    
    axins = inset_axes(ax, width="100%", height="100%", \
                       bbox_to_anchor=((loc[2]-lims[2])/(lims[3]-lims[2]),\
                       (loc[0]-loc[1]/2.0 -lims[0])/(lims[1]-lims[0]), \
                       (loc[3])/(lims[3]-lims[2]), loc[1]/(lims[1]-lims[0])), \
                       bbox_transform=ax.transAxes, borderpad=0)
    axins.set_ylim(0,1)
    axins.set_ylabel('p')
    axins.set_xlim(loc[2],loc[2]+loc[3])
    axins.patch.set_alpha(0.4)
    axins.plot(xs,ps,'-',color=pcolor,linewidth=2)
    axins.get_xaxis().set_visible(False) 
    
axlabels=['a)','b)','c)','d)']

def draw_lines(ax,color2,full=True):
    ax.plot([1900,2100],[.5,.5],'k-.')
    ax.plot([1900,2100],[.159,.159],'k:')
    ax.plot([1900,2100],[.159,.159],color = color2,lw=3,zorder=-1)
    if(full):
        ax.plot([1900,2100],[.841,.841],'k:')
        ax.plot([1900,2100],[.841,.841],color = color2,lw=3,zorder=-1)
    else:
        ax.plot([2028,2100],[.841,.841],'k:')
        ax.plot([2028,2100],[.841,.841],color = color2,lw=3,zorder=-1)
    ax.set_ylim(0,1)
    ax.set_ylabel('P(year)')
    ax.set_xlabel('Year')

fineyears=np.linspace(1950, 2025, 750)
def hatch_below(ax3,thisdates,pvals,color1,shiftup=0):
    if(color1=='indigo'):
        hatchdir="/"
    else:
        hatchdir="\\"
    fineyears=np.linspace(1950, 2100, 1500)
    interpEBMK = np.interp(fineyears,thisdates,pvals) #thres1_ps[:,0]
    sind=np.argmax(interpEBMK>0.159)
    ax3.text(fineyears[sind]-0.1, shiftup, '(',color=color1,ha='center',size=16)
    eind=len(fineyears)- np.argmax(np.flip(interpEBMK)<0.841)
    if (eind!=len(fineyears)):
        ax3.text(fineyears[eind], shiftup, ')',color=color1,ha='center',size=16)
        ax3.plot([fineyears[sind],fineyears[eind]],[shiftup+0.02,shiftup+0.02],color=color1,lw=3)
#        eind=np.argmax(interpEBMK>np.max(interpEBMK)-0.0001)
    xinds=[0]
    xind=np.argmax(interpEBMK>0.5)
    direc=1
    while (xind>0):
        xind=xind+xinds[-1]
        xinds.append(xind)
        direc=-direc
        if(direc==1):
            xind=np.argmax(interpEBMK[xind:]>0.5)
        elif(direc==-1):
            xind=np.argmax(interpEBMK[xind:]<0.5)


    for x in xinds:
        ax3.plot([fineyears[x],fineyears[x]],[shiftup+0.02,interpEBMK[x]],"-",color=color1,alpha=0.5)
        print(fineyears[x])

    #ax3.fill_between(fineyears[sind:eind], np.zeros((eind-sind)), interpEBMK[sind:eind], color=color2,hatch=hatchdir,edgecolor=color1,alpha=0.2)

xh1s=ekf.xh1s
xh0s=ekf.xh0s
stdS=ekf.stdS
stdP=ekf.stdP

#fig, (ax1,ax4)= plt.subplots(2, 1, figsize=(7,10))
fig21, ax1= plt.subplots(1, 1, figsize=(7,6))
fig22, ax4= plt.subplots(1, 1, figsize=(7,6))
#plt.subplots_adjust(top=0.95, bottom=0.05)
ax=ax4
plt.sca(ax)
#plt.plot(dates,xhatminus,'b-',label='a priori EKF estimate', linewidth=3.0)
ax.set_title('Temperature Forecast Probabilities of Threshold Crossings with EBM - Kalman Filter')
xh0s[0]=xh1s[0]
ax.plot(ekf.dates,xh0s,'-.',label='a priori EBM-KF GMST state estimate $\^{T }_{n|n-1}}$', color=ekf.colorekf)

ax.fill_between(ekf.dates, xh0s-2*stdS, xh0s+2*stdS,label="95% CI from $S_n$ of measurements around pred. GMST $\^{T }_{n|n-1}}$", color=ekf.coloruncert)
ax.plot(ekf.dates,ekf.temps,'o',markersize=2,color=ekf.colorgrey,label='noisy HadCRUT5 GMST measurements $Y_{n}$')
lims=[286.1,288.3,1850,2025] #ymin ymax xmin xmax

ekf.plot_boilerplate(ax)
ax.set_xlim([lims[2],lims[3]])
ax.set_ylim([lims[0],lims[1]])
ax.set_yticks(np.arange(286.2,288.3,0.2))
ax.plot([lims[2],lims[3]], [ekf.pindavg, ekf.pindavg], '-',color='sienna',label='thresholds in 0.5K increments')
ax.plot([lims[3]+100,lims[3]+100], [ekf.pindavg, ekf.pindavg], '-',color='purple',linewidth=2,label='P(threshold crossed by given year)')
handles, labels = ax.get_legend_handles_labels()
order = [4,0,1,2,3]
ax.legend([handles[idx] for idx in order],[labels[idx] for idx in order],prop={'size': 9.5})
axb = ax.twinx()
axb.set_yticks(np.arange(12.4,18.4,0.2))
axb.set_ylim(lims[0]- 273.15, lims[1]- 273.15)
axb.set_ylabel('Temperature (°C)')

#plt.tick_params( direction = 'in',bottom=True, top=True, left=True, right=True )#

thres1=ekf.pindavg+0.5
thres2=ekf.pindavg+1.0
thres3=ekf.pindavg+1.5
thres4=ekf.pindavg+2.0
thres5=ekf.pindavg+2.5
t1s=1975
t1l=30
t2s=1998
t2l=30
t3s=2015
t3l=30
t4s=2030
t4l=30
t5s=2045
t5l=30
boxsize=0.5

thres1_ps=np.zeros((ekf.n_iters,1))
thres2_ps=np.zeros((ekf.n_iters,1))
thres1_cps=np.zeros((ekf.n_iters,1))
thres2_cps=np.zeros((ekf.n_iters,1))
for i in range(ekf.n_iters):
    thres1_ps[i]=1-ss.norm.cdf(thres1,xh0s[i],stdS[i])
    thres2_ps[i]=1-ss.norm.cdf(thres2,xh0s[i],stdS[i])
#print(thres1_ps)
#print(thres2_ps)
cdf_subplot(ax,[thres1, boxsize, t1s,t1l],lims,ekf.dates,thres1_ps)
cdf_subplot(ax,[thres2, boxsize, t2s,2025-t2s],lims,ekf.dates,thres2_ps)


#fig, (ax12, ax2) = plt.subplots(2, 1, figsize=(7,10))
fig11, ax12 = plt.subplots(1, 1, figsize=(7,6))
fig12, ax2 = plt.subplots(1, 1, figsize=(7,6))
plt.subplots_adjust(top=0.95, bottom=0.05)
import summary_of_sims as sum_sims
sum_sims.plot_many_sims()
ax2.set_ylim(286.1,288.3)
ax2.set_xlim(1850,2025)
ax2.set_yticks(np.arange(286.2,288.3,0.2))
ax2.set_title('Temperature Forecast Probabilities of Threshold Crossings with CESM2 LENS Model')
ax2.plot([lims[2],lims[3]], [ekf.pindavg, ekf.pindavg], '-',color='sienna',label='thresholds in 0.5K increments')
ax2.plot([lims[3]+100,lims[3]+100], [ekf.pindavg, ekf.pindavg], '-',color='fuchsia',linewidth=2,label='P(threshold crossed by given year)')
ax2.set_ylabel('Temperature (K)')
handles, labels = ax2.get_legend_handles_labels()
order = [3,0,1,2]
ax2.legend([handles[idx] for idx in order],[labels[idx] for idx in order],prop={'size': 9.5})
axb = ax2.twinx()
axb.set_yticks(np.arange(12.4,18.4,0.2))
axb.set_ylim(lims[0]- 273.15, lims[1]- 273.15)
axb.set_ylabel('Temperature (°C)')

#compute moving averages
N = 21
allsimsavgd=np.empty(np.shape(sum_sims.twTR))
for s in range(sum_sims.sms):
    tempssim = sum_sims.twTR[s,:]
    cumsum, moving_aves = [0], []
    for i, x in enumerate(tempssim, 1):
        cumsum.append(cumsum[i-1] + x)
        if i<N/2:
            moving_aves.append(np.nan)
        if i>=N:
            moving_ave = (cumsum[i] - cumsum[i-N])/N
            #can do stuff with moving_ave here
            if(np.isnan(tempssim[i-N:i]).sum()>0):
                moving_ave=np.nan
            moving_aves.append(moving_ave)
    for i in range(int(np.floor(N/2))):
        moving_aves.append(np.nan)
    allsimsavgd[s,:]=moving_aves

ax12.set_title('Climate State Probabilities of Threshold Crossings with CESM2 LENS Model')
#plt.plot(ekf.dates,xh1s,'-',label='a posteriori EKF state estimate $\^{T }_{n}}$', color=ekf.colorekf)
ax12.plot(sum_sims.dates[14:-15],allsimsavgd[0,14:-15],'-',label='single (unforced) 21-year means $_{21}(Y_n)_j$', color=ekf.colorstate, linewidth=0.5,zorder=1)
for r in range(1,sum_sims.sms):
    ax12.plot(sum_sims.dates[14:-15],allsimsavgd[r,14:-15],'-', color=sum_sims.colorgen(r,(70,30,30),(104,195,121)), linewidth=0.5,zorder=1) #(52./255, 235./255, 235./255)
ax12.plot(sum_sims.dates[14:-15],np.nanmean(allsimsavgd,axis=0)[14:-15],'-',label=str(sum_sims.sms)+'-member ensemble average $\overline{_{21}(Y_n)_j}$', color=ekf.colorekf,zorder=3)
ekf.plot_boilerplate(ax12)


ekf.plot_boilerplate(ax12)
ax12.set_xlim([lims[2],lims[3]])
ax12.set_ylim([lims[0],lims[1]])
ax12.set_yticks(np.arange(286.2,288.3,0.2))
axb = ax12.twinx()
axb.set_yticks(np.arange(12.,18.4,0.2))
axb.set_ylim(lims[0]- 273.15, lims[1]- 273.15)
axb.set_ylabel('Temperature (°C)')
ax12.plot([lims[2],lims[3]], [ekf.pindavg, ekf.pindavg], '-',color='sienna',label='thresholds in 0.5K increments')
ax12.plot([lims[3]+100,lims[3]+100], [ekf.pindavg, ekf.pindavg], '-',color='fuchsia',linewidth=2,label='P(threshold crossed by given year)')
handles, labels = ax12.get_legend_handles_labels()
order = [3,0,1,2]
ax12.legend([handles[idx] for idx in order],[labels[idx] for idx in order],prop={'size': 9.5})




thres1_ss=np.zeros((sum_sims.smyrs,1))
thres2_ss=np.zeros((sum_sims.smyrs,1))
thres1_css=np.zeros((sum_sims.smyrs,1))
thres2_css=np.zeros((sum_sims.smyrs,1))
thres3_ss=np.zeros((sum_sims.smyrs,1))
thres4_ss=np.zeros((sum_sims.smyrs,1))
thres5_ss=np.zeros((sum_sims.smyrs,1))
thres3_css=np.zeros((sum_sims.smyrs,1))
thres4_css=np.zeros((sum_sims.smyrs,1))
thres5_css=np.zeros((sum_sims.smyrs,1))
ssn=np.zeros((sum_sims.smyrs,1))

#compute moving avg kalman

Nav = 30
moving_aves_data=np.empty(len(ekf.temps))
moving_aves_data[:] = np.nan
cN=int(np.ceil(Nav/2));
fN=int(np.floor(Nav/2));
for i in range(cN,(len(ekf.temps)-fN+1)):
    lasta=i+fN+1;
    firsta=i-cN;
    moving_aves_data[i] = np.mean(ekf.temps[firsta:lasta]);
fineyears=np.linspace(1950, 2025, 750)
moving_aves_data_interp = np.interp(fineyears,ekf.dates,moving_aves_data)
x1crossavgyr=np.argmax(moving_aves_data_interp>thres1)
#x2crossavgyr=np.argmax(moving_aves_data_interp>thres2)
#print(xcrossavgyr)
plt.figure()
plt.plot(ekf.dates,moving_aves_data)
plt.plot([1900,2000],[thres1,thres1])

plt.legend(prop={'size': 8})
for j in range(sum_sims.smyrs):
    num_sims_now=len([i for i in sum_sims.twTR[:,j] if i > 0])
    thres1_ss[j]=len([i for i in sum_sims.twTR[:,j] if i > thres1])/float(num_sims_now)
    thres2_ss[j]=len([i for i in sum_sims.twTR[:,j] if i > thres2])/float(num_sims_now)
    thres3_ss[j]=len([i for i in sum_sims.twTR[:,j] if i > thres3])/float(num_sims_now)
    thres4_ss[j]=len([i for i in sum_sims.twTR[:,j] if i > thres4])/float(num_sims_now)
    thres5_ss[j]=len([i for i in sum_sims.twTR[:,j] if i > thres5])/float(num_sims_now)
    num_avs_now=len([i for i in allsimsavgd[:,j] if i > 0])
    if (num_avs_now==0):
        num_avs_now=np.nan
    thres1_css[j]=len([i for i in allsimsavgd[:,j] if i > thres1])/float(num_avs_now)
    thres2_css[j]=len([i for i in allsimsavgd[:,j] if i > thres2])/float(num_avs_now)
    thres3_css[j]=len([i for i in allsimsavgd[:,j] if i > thres3])/float(num_avs_now)
    thres4_css[j]=len([i for i in allsimsavgd[:,j] if i > thres4])/float(num_avs_now)
    thres5_css[j]=len([i for i in allsimsavgd[:,j] if i > thres5])/float(num_avs_now)
cdf_subplot(ax2,[thres1, boxsize, t1s,t1l],lims,sum_sims.dates,thres1_ss,'fuchsia')
cdf_subplot(ax2,[thres2, boxsize, t2s,2025-t2s],lims,sum_sims.dates,thres2_ss,'fuchsia')



cdf_subplot(ax12,[thres1, boxsize, t1s,t1l],lims,sum_sims.dates,thres1_css,'fuchsia')
cdf_subplot(ax12,[thres2, boxsize, t2s,2025-t2s],lims,sum_sims.dates,thres2_css,'fuchsia')


thres1_cps=np.zeros((ekf.n_iters,1))
thres2_cps=np.zeros((ekf.n_iters,1))
startf=2022 #not 2023!
endf=2100

for i in range(ekf.n_iters):
    thres1_cps[i]=1-ss.norm.cdf(thres1,xh1s[i],stdP[i])
    thres2_cps[i]=1-ss.norm.cdf(thres2,xh1s[i],stdP[i])

plt.rcParams['figure.figsize'] = (7,7)



fdates=np.arange(startf,endf,1/4)
fdata = np.load("SSP370.npz")
dif_array=np.abs(fdata['gmstheights']-thres3)
incp3=dif_array.argmin()
dif_array=np.abs(fdata['gmstheights']-thres4)
incp4=dif_array.argmin()
dif_array=np.abs(fdata['gmstheights']-thres5)
incp5=dif_array.argmin()
fcss=(1000-fdata['cgmst'])/1000
ftpss=(1000-fdata['tpgmst'])/1000
sup=0.05

fig2, axes=plt.subplots(nrows=2, ncols=5,figsize=(14,8))
fig2.subplots_adjust(wspace=0.45,hspace=0.45)

if(True):
    data3 = np.genfromtxt(open("KF6projectionSSP.csv", "rb"),dtype=float, delimiter=',')
    (fxhat,fP0)=ekf.ekf_future(startf,ekf.Plastretain,ekf.xlastretain,endf+1,np.log10(data3[:,1+2]),data3[:,6+2],noVolcs=True)
    fxh1s=fxhat[:,0]
    fstdP=np.sqrt(np.abs(fP0))[:,0,0]
    fstdS=np.sqrt(np.abs(fP0+ekf.Rtvara))[:,0,0]
    thres3_fcps=np.zeros((endf-startf,1))
    thres4_fcps=np.zeros((endf-startf,1))
    thres5_fcps=np.zeros((endf-startf,1))
    thres3_fps=np.zeros((endf-startf,1))
    thres4_fps=np.zeros((endf-startf,1))
    thres5_fps=np.zeros((endf-startf,1))
    for i in range(endf-startf):
        thres3_fcps[i]=1-ss.norm.cdf(thres3,fxh1s[i],fstdP[i])
        thres4_fcps[i]=1-ss.norm.cdf(thres4,fxh1s[i],fstdP[i])
        thres5_fcps[i]=1-ss.norm.cdf(thres5,fxh1s[i],fstdP[i])
        thres3_fps[i]=1-ss.norm.cdf(thres3,fxh1s[i],fstdS[i])
        thres4_fps[i]=1-ss.norm.cdf(thres4,fxh1s[i],fstdS[i])
        thres5_fps[i]=1-ss.norm.cdf(thres5,fxh1s[i],fstdS[i])
    fdate0=np.arange(startf,endf,1)
    axes[0,4].plot(fdate0,thres5_fcps,'-',color='blue',linewidth=2)
    hatch_below(axes[0,4],fdate0,thres5_fcps[:,0],'blue',sup/2)
    axes[1,4].plot(fdate0,thres5_fps,'-',color='blue',linewidth=2)
    hatch_below(axes[1,4],fdate0,thres5_fps[:,0],'blue',sup/2)
    axes[0,3].plot(fdate0,thres4_fcps,'-',color='blue',linewidth=2)
    hatch_below(axes[0,3],fdate0,thres4_fcps[:,0],'blue',sup/2)
    axes[1,3].plot(fdate0,thres4_fps,'-',color='blue',linewidth=2)
    hatch_below(axes[1,3],fdate0,thres4_fps[:,0],'blue',sup/2)
    axes[0,2].plot(fdate0,thres3_fcps,'-',color='blue',linewidth=2)
    hatch_below(axes[0,2],fdate0,thres3_fcps[:,0],'blue',sup/2)
    axes[1,2].plot(fdate0,thres3_fps,'-',color='blue',linewidth=2,label="EBM-KF Unif")
    hatch_below(axes[1,2],fdate0,thres3_fps[:,0],'blue',sup/2)
    
    

ax3 = axes[0,4]
draw_lines(ax3,ekf.colorstate)
ax3.set_xlim(t5s,t5s+t4l)
ax3.plot(fdates,fcss[incp5,:],'-',color='indigo',linewidth=2, label='Extended Kalman Filter on Measurements')
ax3.plot(sum_sims.dates,thres5_css,'-',color='fuchsia',linewidth=2, label='CESM2 LENS Model Simulations')
ax3.set_title('e)          Climate State          \n +2.5K (SSP370)')
hatch_below(ax3,fdates,fcss[incp5,:],'indigo')
hatch_below(ax3,sum_sims.dates,thres5_css[:,0],'fuchsia',sup)
minor_xticks = np.arange(t5s,t5s+t5l, 5)
ax3.set_xticks(minor_xticks, minor=True)

ax3 = axes[1,4]
draw_lines(ax3,ekf.coloruncert,)
ax3.set_xlim(t5s,t5s+t4l)
ax3.plot(fdates,ftpss[incp5,:],'-',color='indigo',linewidth=2, label='EBM - Kalman')
ax3.plot(sum_sims.dates,thres5_ss,'-',color='fuchsia',linewidth=2, label='LENS2 Sims')
ax3.set_title('j)    Temperature Forecast     \n +2.5K (SSP370)')
hatch_below(ax3,fdates,ftpss[incp5,:],'indigo')
hatch_below(ax3,sum_sims.dates,thres5_ss[:,0],'fuchsia',sup)
minor_xticks = np.arange(t5s,t5s+t5l, 5)
ax3.set_xticks(minor_xticks, minor=True)


ax3 = axes[0,3]
draw_lines(ax3,ekf.colorstate)
ax3.set_xlim(t4s,t4s+t4l)
ax3.plot(fdates,fcss[incp4,:],'-',color='indigo',linewidth=2, label='Extended Kalman Filter on Measurements')
ax3.plot(sum_sims.dates,thres4_css,'-',color='fuchsia',linewidth=2, label='CESM2 LENS Model Simulations')
ax3.set_title('d)          Climate State          \n +2.0K (SSP370)')
hatch_below(ax3,fdates,fcss[incp4,:],'indigo')
hatch_below(ax3,sum_sims.dates,thres4_css[:,0],'fuchsia',sup)
minor_xticks = np.arange(t4s,t4s+t4l, 5)
ax3.set_xticks(minor_xticks, minor=True)

ax3 = axes[1,3]
draw_lines(ax3,ekf.coloruncert)
ax3.set_xlim(t4s,t4s+t4l)
ax3.plot(fdates,ftpss[incp4,:],'-',color='indigo',linewidth=2, label='EBM - Kalman')
ax3.plot(sum_sims.dates,thres4_ss,'-',color='fuchsia',linewidth=2, label='LENS2 Sims')
ax3.set_title('i)    Temperature Forecast     \n +2.0K (SSP370)')
hatch_below(ax3,fdates,ftpss[incp4,:],'indigo')
hatch_below(ax3,sum_sims.dates,thres4_ss[:,0],'fuchsia',sup)
minor_xticks = np.arange(t4s,t4s+t4l, 5)
ax3.set_xticks(minor_xticks, minor=True)

ax3 = axes[0,2]
draw_lines(ax3,ekf.colorstate)
ax3.set_xlim(t3s,t3s+t3l)
ax3.plot(fdates,fcss[incp3,:],'-',color='indigo',linewidth=2, label='Extended Kalman Filter on Measurements')
ax3.plot(sum_sims.dates,thres3_css,'-',color='fuchsia',linewidth=2, label='CESM2 LENS Model Simulations')
ax3.set_title('c)          Climate State          \n +1.5K (SSP370)')
hatch_below(ax3,fdates,fcss[incp3,:],'indigo')
hatch_below(ax3,sum_sims.dates,thres3_css[:,0],'fuchsia',sup)
minor_xticks = np.arange(t3s,t3s+t3l, 5)
ax3.set_xticks(minor_xticks, minor=True)

ax3 = axes[1,2]
draw_lines(ax3,ekf.coloruncert,full=False)
ax3.set_xlim(t3s,t3s+t3l)
ax3.plot(fdates,ftpss[incp3,:],'-',color='indigo',linewidth=2, label='EBM-KF Volc')
ax3.plot(sum_sims.dates,thres3_ss,'-',color='fuchsia',linewidth=2, label='LENS2')
ax3.set_title('h)    Temperature Forecast     \n +1.5K (SSP370)')
hatch_below(ax3,fdates,ftpss[incp3,:],'indigo')
hatch_below(ax3,sum_sims.dates,thres3_ss[:,0],'fuchsia', sup)
minor_xticks = np.arange(t3s,t3s+t3l, 5)
ax3.set_xticks(minor_xticks, minor=True)
ax3.legend(loc='upper right', prop={'size': 8})
ax3.legend_.set_bbox_to_anchor((0.44, 1.00))

######
ax3 = axes[0,1]
draw_lines(ax3,ekf.colorstate)
ax3.set_xlim(t2s,t2s+t2l)
ax3.plot(ekf.dates,thres2_cps,'-',color='indigo',linewidth=2, label='Extended Kalman Filter on Measurements')
ax3.plot(sum_sims.dates,thres2_css,'-',color='fuchsia',linewidth=2, label='CESM2 LENS Model Simulations')
ax3.set_title('b)          Climate State          \n +1.0K Threshold')
hatch_below(ax3,ekf.dates,thres2_cps[:,0],'indigo')
hatch_below(ax3,sum_sims.dates,thres2_css[:,0],'fuchsia',sup)
minor_xticks = np.arange(1995,t2s+t2l, 5)
ax3.set_xticks(minor_xticks, minor=True)

ax3 = axes[1,1]
draw_lines(ax3,ekf.coloruncert)
ax3.set_xlim(t2s,t2s+t2l)
ax3.plot(ekf.dates,thres2_ps,'-',color='indigo',linewidth=2, label='EBM - Kalman')
ax3.plot(sum_sims.dates,thres2_ss,'-',color='fuchsia',linewidth=2, label='LENS2 Sims')
ax3.set_title('g)    Temperature Forecast     \n +1.0K Threshold')
hatch_below(ax3,ekf.dates,thres2_ps[:,0],'indigo')
hatch_below(ax3,sum_sims.dates,thres2_ss[:,0],'fuchsia',sup)
#ax3.text(fineyears[x2crossavgyr],0,'$\\ast$',color='k',ha='center',size=14)
minor_xticks = np.arange(1995,t2s+t2l, 5)
ax3.set_xticks(minor_xticks, minor=True)

ax3 = axes[0,0]
draw_lines(ax3,ekf.colorstate)
ax3.set_xlim(t1s,t1s+t1l)
ax3.plot(ekf.dates,thres1_cps,'-',color='indigo',linewidth=2, label='Extended Kalman Filter on Measurements')
ax3.plot(sum_sims.dates,thres1_css,'-',color='fuchsia',linewidth=2, label='CESM2 LENS Model Simulations')
ax3.set_title('a)          Climate State          \n +0.5K Threshold')
hatch_below(ax3,ekf.dates,thres1_cps[:,0],'indigo')
hatch_below(ax3,sum_sims.dates,thres1_css[:,0],'fuchsia',sup)
minor_xticks = np.arange(t1s,t1s+t1l, 5)
ax3.set_xticks(minor_xticks, minor=True)

ax3.text(fineyears[x1crossavgyr],0.5,'$\\ast$',color='darkorange',ha='center',size=20)

ax3 = axes[1,0]
draw_lines(ax3,ekf.coloruncert)
ax3.set_xlim(t1s,t1s+t1l)
ax3.plot(ekf.dates,thres1_ps,'-',color='indigo',linewidth=2, label='EBM-KF')
ax3.plot(sum_sims.dates,thres1_ss,'-',color='fuchsia',linewidth=2, label='LENS2')
ax3.set_title('f)    Temperature Forecast     \n +0.5K Threshold')
hatch_below(ax3,ekf.dates,thres1_ps[:,0],'indigo')
hatch_below(ax3,sum_sims.dates,thres1_ss[:,0],'fuchsia', sup)
minor_xticks = np.arange(t1s,t1s+t1l, 5)
ax3.set_xticks(minor_xticks, minor=True)
ax3.legend(loc='best', prop={'size': 8})
ax3.text(fineyears[x1crossavgyr],0.5,'$\\ast$',color='darkorange',ha='center',size=20)






ax=ax1
plt.sca(ax)
#plt.plot(dates,xhatminus,'b-',label='a priori EKF estimate', linewidth=3.0)
plt.title('Climate State Probabilities of Threshold Crossings with EBM - Kalman Filter')
plt.plot(ekf.dates,xh1s,'-',label='a posteriori EBM-KF GMST state estimate $\^{T }_{n}}$', color=ekf.colorekf)

#plt.fill_between(ekf.dates, xh1s-stdP, xh1s+stdP,label="68% CI of EKF state ($\sqrt{ P_{n}}$)", color=ekf.colorstate)
plt.fill_between(ekf.dates, xh1s-2*stdP, xh1s+2*stdP,label="95% CI from $P_n$ of GMST state $\^{T }_{n}}$", color=ekf.colorstate,alpha=0.3,zorder=0)
plt.plot(ekf.dates,ekf.temps,'o',markersize=2,color=ekf.colorgrey,label='noisy HadCRUT5 GMST measurements $Y_{n}$')

ekf.plot_boilerplate(ax)
ax.set_xlim([lims[2],lims[3]])
ax.set_ylim([lims[0],lims[1]])
ax.set_yticks(np.arange(286.2,288.3,0.2))

ax.plot([lims[2],lims[3]], [ekf.pindavg, ekf.pindavg], '-',color='sienna',label='thresholds in 0.5K increments')
ax.plot([lims[3]+100,lims[3]+100], [ekf.pindavg, ekf.pindavg], '-',color='purple',linewidth=2,label='P(threshold crossed by given year)')
handles, labels = plt.gca().get_legend_handles_labels()
order = [4,0,1,2,3]
plt.legend([handles[idx] for idx in order],[labels[idx] for idx in order],prop={'size': 9.5})

plt.tick_params( direction = 'in',bottom=True, top=True, left=True, right=True )
axb = ax.twinx()
axb.set_yticks(np.arange(12.4,18.4,0.2))
axb.set_ylim(lims[0]- 273.15, lims[1]- 273.15)
axb.set_ylabel('Temperature (°C)')
#print(thres1_ps)
#print(thres2_ps)

cdf_subplot(ax,[thres1, boxsize, t1s,t1l],lims,ekf.dates,thres1_cps)
cdf_subplot(ax,[thres2, boxsize, t2s,2025-t2s],lims,ekf.dates,thres2_cps)


fig2.savefig("cdf4.pdf",format="pdf")
fig2.savefig("cdf4.png", dpi=400,format="png")
fig11.savefig("cdf1_lens.pdf",format="pdf")
fig11.savefig("cdf1_lens.png", dpi=400,format="png")
fig12.savefig("cdf2_lens.pdf",format="pdf")
fig12.savefig("cdf2_lens.png", dpi=400,format="png")
fig21.savefig("cdf1_ekf.pdf",format="pdf")
fig21.savefig("cdf1_ekf.png", dpi=400,format="png")
fig22.savefig("cdf2_ekf.pdf",format="pdf")
fig22.savefig("cdf2_ekf.png", dpi=400,format="png")

plt.show()
